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1  Introduction 


The  focus  of  this  AFOSR  sponsored  research  project  is  to  investigate  the  small-scale  (<  100  nr) 
plasma  dynamics  of  Convective  Equatorial  Ionospheric  Storms  (CEIS),  i.e.  Equatorial  Spread- 
F.  CEIS  are  challenging  to  predict,  and  they  cause  debilitating  degradation  of  earth-to-satellite 
communication  signals  (e.g.  GPS).  The  instabilities  at  the  smallest  length  scale  are  essential  for 
rationalizing  the  existence  of  coherent  radar  echos  of  CEIS,  and  these  instabilities  provide  the 
energetic  sink  for  the  large  scale  dynamics,  commonly  called  spread-F  bubbles.  Little  direct  evidence 
exists  to  identify  the  specific  driving  force  responsible  for  developing  these  waves,  and  several 
proposed  theories  for  their  growth  mechanism  fail  to  address  all  length  scales  [Huba  and  Ossakow, 
1979,  1981a, b].  Our  objective  was  to  develop  a  better  understanding  of  this  physics  by  simulating 
it  using  Particle-in-Cell  (PIC)  codes. 

PIC  simulators  model  small-scale,  non-linear,  and  kinetic  effects,  which  the  large-scale  fluid- 
based  simulations  cannot  study.  However,  PIC  codes  cannot  explore  large-scale  or  slowly  evolving 
systems.  Nevertheless,  we  used  such  a  code  to  investigate  2  and  3  dimensional  systems  initialized 
with  strong  density  gradients.  We  found  that  these  gradients  resulted  in  small  scale  drift  wave 
instabilities  perpendicular  to  the  background  density  gradient.  The  gradient  length  scales  required 
to  produce  these  are  slightly  smaller  than  those  previously  measured  in  situ  by  satellites  and  rockets. 
Therefore,  we  conclude  that  a  simple  drift  wave  instability  is  sufficient  to  explain  the  development 
of  small  scale  irregularities  so  long  as  they  quickly  dissipate  the  steepened  structures  that  drive 
their  growth.  We  found  that  the  third  dimension  (dynamics  parallel  to  the  background  magnetic 
field,  B)  resulted  in  subtly  different  dynamics.  Since  the  dynamics  developed  much  slower  parallel 
to  B,  the  2  dimensional  simulations,  perpendicular  to  B,  produced  qualitatively  similar  short-time 
dynamics  and  the  2D  simulations  are  still  very  relevant.  We  also  performed  exploratory  studies 
into  the  combined  driver  of  local  electric  fields  and  density  gradients.  We  are  currently  preparing 
a  manuscript  discussing  these  results. 

While  these  PIC  simulations  are  well  suited  to  accurately  describe  the  physics  at  these  small 
scales,  the  calculations  are  computationally  costly  compared  to  fluid-based  models.  Considerable 
effort  for  this  project  was  directed  towards  improving  the  capabilities  of  our  PIC  simulator  to  better 
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suit  it  for  these  large  simulations.  The  major  update  involved  adding  the  capability  non-periodic 
simulations,  where  the  flux  of  particles  at  one  simulation  boundary  differed  from  the  opposite 
boundary.  This  ability  made  it  possible  to  sustain,  essentially  indefinitely,  density  gradients  across 
the  simulation.  The  fact  that  this  new  code  is  parallelized  is  also  important,  because  it  allows  us 
to  harness  the  efficiency  of  some  of  the  worlds  largest  super  computers  and  to  simulate  for  larger 
length  scales  and  longer  times.  We  also  implemented  a  quasi-neutral  solver  that  treats  the  electrons 
as  an  inertialless  fluid  which  responds  instantly  to  the  ion  density  while  maintaining  V  •  J  =  0. 
While  this  formulation  has  less  sever  numerical  constraints,  it  is  much  more  difficult  to  parallelize. 

2  Detailed  Discussion 

Density  gradients  play  an  essential  role  in  driving  CEIS,  and  they  are  the  main  source  of  energy 
for  the  largest  scale  irregularities.  As  in  situ  measurements  have  verified,  steep  local  density  gra¬ 
dients  exist  following  the  development  of  spread-F  bubbles.  The  spectra  of  density  irregularities 
perpendicular  to  the  magnetic  field  differ  horizontally  vs.  vertically,  and  particularly  sharp  gra¬ 
dients  (gradient  length  scales-n(Vn)-1,  approaching  50  m)  have  been  measured  vertically  [Singh 
and  Szuszczewicz,  1984],  The  theory  of  collisional  interchange  instability  explains  this  steepening 
preferentially  along  the  vertical  direction  [Kelley  et  al.,  1987].  Hysell  et  al.  [1994]  showed  that  for 
a  simple  model  of  these  waves,  parallel  to  the  background  density  gradient  (typically  the  vertical 
direction)  they  lead  to  steepened  intermediate  scale  (100-1000  m)  structures  with  very  sharp  gra¬ 
dients.  It  is  clear  from  measurement  and  theory  that  these  steep  intermediate  scale  gradients  exist, 
and  it  is  these  structures  (and  their  associated  fields)  we  hypothesize  predominately  drive  the  small 
scale  irregularities. 

2D  simulations  Our  2D  simulations  were  oriented  perpendicular  to  the  magnetic  field  and  en¬ 
compassed  approximately  2500  nr2.  Since  the  time  to  develop  the  intermediate  scale  structures 
(nearly  15  minutes)  is  much  longer  then  can  be  physically  modeled  using  PIC  simulations,  we  ini¬ 
tialize  the  system  with  a  density  gradient  comparable  to  that  which  would  be  produced  by  another 
process,  such  as  the  collisional  interchange  instability  [Hysell  et  al.,  1994],  In  the  direction  perpen- 
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Parameters 

Value 

Sx,  5y 

nr 

0.06 

Nx 

1024 

Ny 

512 

St 

s 

2.4e-06 

B0  ■  z 

G 

2.500e-05 

nproc 

256 

L  =  n(  Vn  )_1 

nr 

7.680e+00 

n0 

m 

3.078e+08 

To 

K 

1.000e+03 

me 

kg 

9.100e-29 

N particle ,  e 

5.244e+07 

^ n 

s-1 

0.000e+00 

Cyclotron  frequency,  fie 

s-1 

4.401e+04 

Larmour  radius,  ae 

m 

2.799e-01 

m,; 

kg 

2.600e-26 

N particle ,  i 

1.049e+08 

no 

m 

3.078e+08 

Cyclotron  frequency,  11, 

s-1 

1.540e+02 

Larmour  radius,  a. 

m 

4.731e+00 

Table  1:  Simulation  Parameters 


dicular  to  the  gradient  (y),  we  used  periodic  boundary  conditions  while  parallel  (x)  we  use  more  a 
range  of  non-periodic  boundary  conditions,  discussed  below.  Table  1  shows  the  parameters  used  to 
initialize  simulator.  We  used  values  typical  of  nighttime  F-region  conditions,  with  two  exceptions: 
The  peak  density  is  the  lowest  value  typically  found  during  the  nighttime.  With  this  value,  the 
plasma  Debye  length  is  largest,  resulting  in  the  most  efficient  simulation  possible.  Comparison 
runs  at  larger  densities  produced  similar  results.  Secondly,  the  electron  mass  is  larger  than  the 
physical  value,  and  this  too  is  chosen  to  maximize  computational  resources.  Using  heavy  electrons 
is  a  common  approach  in  PIC  simulations,  and  we  stayed  within  suggested  limits  [Dirnant  and 
Oppenheim,  2004],  The  electron’s  mass’s  effect  was  studied  and  is  discussed  below. 

In  such  a  system  one  would  expect  diamagnetic  drifts  to  develop,  where  ions  drift  upward 
and  electrons  downward.  This  drives  a  drift  wave  instability,  which  has  been  studied  extensively 
[Gary,  1993].  This  instability,  and  similar  variants,  are  stable  for  the  weak  gradients  that  drive  the 
large  scale  CEIS  dynamics  [Huba  and  Ossakow,  1981b],  but  the  standard  analysis  assumes  that 
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the  gradient  length  scale  (n(Vn)  x)  exceeds  the  ion  Larmour  radius.  This  assumption  fails  for  the 
steep  structures  we  used.  Figure  1  shows  a  typical  evolution  of  the  ion  density. 


nx{x,y,t:  dt  =15.0528  ms)fnp€ak 


x  (m) 

Figure  1:  The  ion  density  as  a  function  of  time,  for  a  simulation  with  the  parameters  listed  in  Table 
1.  Time  progress  from  left  to  right  and  top  to  bottom,  and  each  frame  shows  high  density  (red)  on 
the  right  and  low  density  on  the  left,  with  a  transition  in  the  middle. 

Early  on  in  the  simulation,  regular,  sinusoidal  drift  waves  form  as  seen  in  the  third  pane  in 
Figure  1).  As  time  progresses,  the  linear  process  gives  way  to  a  non-linear  turbulent  dynamics, 
still  localized  to  the  region  with  the  highest  gradient.  Surprisingly,  the  turbulence  takes  the  form 
of  eddies.  A  visual  inspection  of  density  perturbations  suggests  that  the  longest  wave  lengths  have 
the  largest  growth  rates,  and  a  spectral  analysis  confirms  that  the  instability  consists  of  modes 
whose  wave  lengths  are  larger  than  1  m. 
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Figure  2:  A  spectral  decomposition  of  the  ion  density  shown  as  ky  vs  time  for  kx  =  0. 


Instability  Growth  and  Saturation  Figure  2  presents  the  amplitude  of  each  ky  (kx  =  0)  mode 
as  a  function  of  time,  confirming  that  waves  smaller  than  1  m  (a  numerical  limitation???)  are 
effectively  stable.  It  also  shows  in  the  left  most  plot  that  1,  3  and  6  nr  waves  grow  simultaneously, 
with  a  similar  growth  rate. 
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.  Ion  Drift  ♦ — •  e-  Diamagnetic  Drift 

--  e- Drift  Ambipolar  ExB  Drift 

—  3  m  Phase  Velocity  *  *  ExB+Diamag.  Drifts 


Figure  3:  Physical  and  theoretical  drifts  as  a  function  of  the  gradient  direction  (x)  and  averaged  in 
the  perpendicular  direction  (y). 


The  diamagnetic  “drift”  results  from  more  particles  gyrating  in  the  higher  density  region,  so 
that  on  average  more  ions  move  vertically  while  on  average  more  electrons  descend.  We  start  our 
particles,  however,  with  zero  drift,  and  through  the  conservation  of  momentum  in  a  collisionless 
plasma,  the  ions  do  not  drift  significantly  in  the  actual  simulations.  On  the  other  hand,  the 
electrons  drift  with  twice  their  expected  velocity.  The  extra  driver  takes  the  form  of  an  anomalous 
ambipolar  field  that  develops  to  counter  the  ions  response  to  a  pressure  driven  Lorentz  force,  and 
thus  the  field  helps  effectively  conserve  momentum.  On  the  other  hand,  the  ambipolar- ExB  drift 
is  coincident  with  the  electrons  diamagnetic  drift,  and  the  electrons  combined  drift  is  twice  as 
fast.  If  we  initialize  the  particles  with  their  corresponding  drifts  at  t=0,  then  the  ambipolar  field 
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does  not  develop,  but  the  particles  still  have  a  relative  drift  equal  to  the  difference  between  their 
diamagnetic  drift  velocities.  The  instability  progresses  with  the  ions  and  at  velocities  smaller  than 
the  ion  acoustic  speed. 


Figure  4:  n(ky,uj),  the  spectra  for  the  ion  distribution. 


By  investigating  the  density  irregularities  in  Fourier  space,  we  can  identify  trends  that  may 
suggest  a  named  instability  and  the  results  at  fixed  wave  length  are  comparable  to  the  coherent 
echoes  measured  by  ionospheric  radars.  Figure  4  presents  such  an  analysis,  showing  the  ion  density 
as  a  function  of  ky  and  u).  Its  interesting  to  see  that  for  small  ky,  these  small  wave  length  modes 
have  insignificant  amplitudes,  which  is  equivalent  to  that  which  is  plotted  in  Figure  2.  Also  of 
interest  is  the  dispersion  curve  for  the  modes  with  wave  lengths  larger  than  1  m  (ky  ~  6.3;  the 
frequency  (w)  has  a  k2  dependence,  similar  in  nature  to  diffusion  (V2n  — >  k2n  ).  Hysell  and  Kelley 
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[1997]  suggested  that  the  small  scale  waves  are  likely  damped  by  diffusion,  and  such  a  process  would 
have  a  power  spectrum  comparable  to  those  measured  by  rockets  and  satellites.  Since  there  are  no 
collisions  in  these  simulations,  classical  diffusion  is  not  possible,  but  the  resulting  dynamics  has  a 
diffusion-like  behavior  that  is  consistent  with  the  analysis  of  Hysell  and  Kelley  [1997]. 

We  investigated  a  range  of  gradients  that  produced  these  instabilities,  and  we  found  that  for 
gradient  length  scales  larger  than  22  nr,  the  instability  did  not  form.  This  is  at  least  half  as  long 
(50  nr)  as  the  steepest  gradients  measured  in  situ.  This  may  suggest  that  these  steep  gradients  are 
transient  in  CEIS,  and  this  instability  provides  the  mechanism  for  rapidly  dissipating  such  steep 
gradients.  The  simulated  instability  quickly  grows  and  saturates  (see  Figure  2)  within  a  fraction  of 
a  second  (.1  to  .5  seconds,  depending  on  the  gradient  strength).  The  turbulent  dynamics  facilitates 
the  dissipation  of  such  density  gradients.  In  CEIS,  the  large  scale  background  gradient  (along 
with  other  forces)  drives  the  growth  of  large  scale  (>  1  knr)  structures  that  develop  progressively 
steeper  gradients  until  this  small  scale  drift  wave  instability  can  rapidly  grow  and  dissipate  the 
medium-scale  (100-1000  nr)  waves. 

Boundary  conditions,  simulation  size  and  electron  mass  This  grant  supported  several 
key  improvements  to  our  simulation  code  that  enabled  us  to  study  much  larger,  non-periodic 
simulations.  We  implemented  a  Message  Passing  Interface  (MPI)  parallel  field  solver,  extending  an 
existing  serial  version  of  the  code  [Birdsall  and  Langdon,  1991,  Austin  et  al.,  2004].  This  allowed 
us  to  efficiently  deploy  our  simulator  with  as  many  as  16384  processors.  This  parallelization  is 
essential  to  try  to  simulate  as  broad  a  range  of  length  scales  as  possible,  and  since  the  lower  end  of 
the  range  (approximately  1  cm)  is  fixed  by  numerical  requirements,  the  added  parallelization  is  the 
only  way  to  reach  the  scales  relevant  for  this  study.  The  non-periodic  conditions  allow  for  several 
possible  field  and  particle-injection  boundary  conditions. 

For  the  field,  Dirichlet  (0  =  0),  Neumann  (E  =  0),  or  Open  (0(oo)  =  0)  boundary  conditions  are 
available,  and  the  particles  are  injected  from  a  fixed  flux  or  reflected  at  the  boundary  to  maintain 
a  fixed  density  and  temperature.  The  system  responds  to  excess  change,  so  it  is  essential  to  use 
the  correct  boundary  conditions,  otherwise  instabilities  tend  to  develop  on  the  boundary.  Our 
experimentation  with  various  combinations  concluded  that  either  Dirichlet  or  Neumann  boundary 


conditions  provided  the  most  stable  representation  of  the  electric  potential,  and  while  both  injection 
and  reflection  particle  boundary  conditions  provide  similar  results,  reflection  boundary  conditions 
are  less  likely  to  result  in  charge  build  up  on  the  boundaries. 

For  our  studies,  we  investigated  a  range  of  gradients,  and  in  comparison  studies,  we  maintained 
all  parameters  except  for  the  initial  density  profile.  The  distance  over  which  the  density  changed 
remained  fixed  as  did  the  peak  density,  so  the  lowest  density  changed  with  the  gradient  strength. 
Could  it  be  that  by  having  a  gradient  length  scale  longer  than  the  distance  over  which  the  density 
changed  affect  the  nature  of  the  instability?  For  example,  our  densities  changed  over  15.5  m  in 
most  simulations,  but  we  studied  gradient  length  scales  as  much  as  30  m.  For  these  conditions,  no 
waves  developed,  so  is  this  due  to  the  limited  range?  When  we  extended  the  range  to  60  nr,  we  still 
found  that  the  waves  were  stable,  from  which  we  are  able  to  confidently  conclude  that  the  cutoff 
for  instability  growth  is  approximately  22  nr.  While  we  would  ideally  wish  to  study  much  larger 
systems,  we  found  our  limited  size  did  not  prevent  us  from  performing  informative  simulations. 

Another  parameter  that  could  influence  our  simulations  is  the  electron  mass,  which  is  artificially 
increased  in  the  simulations  to  decrease  the  severity  of  certain  numerical  requirements.  Upon 
lowering  the  mass  to  value  approaching  the  true  electron  mass,  the  simulation  results  remained 
qualitatively  the  same.  One  key  difference  is  a  subtle  increase  in  the  growth  of  very  short  wavelength 
(<  0.5  nr)  modes,  but  the  longest  wavelength  (1  to  6  nr)  modes  grew  at  the  same  rate  and  saturated 
at  the  same  amplitude.  We  conclude  that  kinetic  effects  dissipate  the  smallest  wavelength  modes  in 
the  larger-nrass  simulations,  and  these  effects  are  dependent  on  the  electron  mass.  The  qualitative 
finding  that  the  largest  waves  grow  due  to  the  gradient  does  not  change.  We  are  still  exploring  which 
kinetic  process  specifically  contributes  to  the  decay  of  the  smallest  modes  in  these  simulations. 

Finally,  we  considered  the  implications  of  using  a  regular  grid,  in  space  and  time,  and  we  were 
conscious  that  frequently  these  numerical  approximations  lead  to  numerical  instabilities  that  are 
challenging  to  distinguish  from  physical  processes.  Not  only  do  all  essential  physical  time  (plasma 
oscillation  period,  electron  gyro  period,  etc)  and  length  (e.g.  Debye  length)  scales  need  to  be 
resolved,  they  must  be  resolved  several  times  over.  Specifically,  the  simulation  parameters  resulted 
in  a  relative  large  electron  cyclotron  frequency  (all  though  smaller  than  the  physical  frequency), 
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and  the  simulations  that  did  not  sufficiently  resolve  the  entire  gyration  ultimately  lead  to  heating. 
Just  resolving  this  cyclotron  frequency  leads  to  a  1000  K  rise  in  the  temperature,  but  this  heating 
was  eliminated  if  the  time  step  became  1  percent  of  the  gyration  period.  This  of  course  made  each 
calculation  100  fold  more  expensive  than  originally  perceived,  but  that  is  the  cost  of  such  reliable 
results. 

We  also  explored  developing  a  parallel  quasi-neutral  hybrid  solver.  Using  the  assumptions 
nt  =  ne,  V  •  J  =  0,  and  inertialless  electrons,  the  electron  momentum  equation  is  solved  for  the 
electric  potential.  Making  numerical  approximations  leads  to  a  non-nornral  matrix  form  that  is  slow 
to  solve,  but  more  importantly,  the  problem  scales  poorly.  In  spite  of  employing  parallel  matrix 
solvers,  the  code  only  scaled  to  10’s  of  processors,  much  less  than  the  pure-PIC  implementation. 

3  dimensional  simulations  For  our  3D  simulations,  we  studied  a  system  with  simulation  box 
dimensions  120  nr  by  60  nr  by  30  nr  for  the  magnetic  field  (B),  density  gradient,  and  perpendicular 
directions  respectively.  Perpendicular  to  B  these  conditions  are  identical  to  those  chosen  for  the 
2D  simulations,  and  all  but  the  gradient  direction  are  treated  with  periodic  boundary  conditions. 
We  found  very  little  density  perturbations  with  wave  vector  components  parallel  to  B,  but  this  may 
have  been  limited  by  two  factors.  First  it  is  expected  that  the  instability  along  B  should  develop 
much  longer  wavelength  modes,  and  the  simulation  box  dimensions  along  B  may  be  too  small 
to  resolve  the  fastest  growing  parallel  modes.  Secondly,  the  development  for  all  modes  parallel  to 
B  is  likely  slower  than  the  perpendicular  direction.  With  this  in  mind,  we  feel  confident  that  the 
2D  simulations,  which  produce  nearly  identical  results,  provide  a  good  representation  of,  at  least, 
the  instabilities  development-time  dynamics. 

2D  bubble  We  also  explored  the  development  of  drift  waves  for  a  circular  density  depletion,  that 
is  a  miniature  spread-F  bubble.  We  employed  periodic  boundary  conditions  for  both  directions, 
but  otherwise  the  plasma  parameters  are  the  same  as  used  above.  In  the  intermediate  region, 
the  gradient  is  constant,  as  opposed  to  the  exponential  transition  (to  impose  a  constant  gradient 
length  scale)  used  in  the  non-periodic  simulations  presented  above.  This  results  in  a  pressure  that 
decreases  away  from  the  depletions  center,  and  therefore  the  instability  grows  predominantly  in  the 
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lower  density  end  of  the  gradient.  Otherwise  the  dynamics  is  similar,  shown  in  a  different  format  in 
Figure  5.  Since  the  pressure  is  stronger  near  the  center,  the  electrons  drift  faster  around  the  bubble, 
and  ultimately  this  faster  drift  induces  the  slowly  propagating  waves  to  travel  slightly  more  quickly 
in  the  center  as  well.  The  difference  produces  a  spiral  pattern,  but  it  also  leads  to  turbulence  sooner 
in  the  center.  The  challenge  with  periodic  simulations  is  that  the  driving  density  gradient  dissipates 
before  the  instability  can  completely  saturate.  This  problem,  however,  reinforces  the  conclusion 
that  the  steep  gradients  used  in  these  simulations  are  not  observed  because  they  are  transient  in 
the  F-region. 


Density  Difference:  —  n*on  (*  =  o) 
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Figure  5:  The  density  difference  (n(t) 


—  n(t  =  0)  for  a  periodic,  circular  density  depletion. 


We  also  briefly  investigated  the  effect  of  an  ambient  electric  field  on  the  instability’s  develop¬ 
ment,  because  satellites,  such  as  C/NOFs,  have  measured  electric  fields  associated  with  the  walls 
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of  the  bubbles.  Figure  6  shows  a  simulation  with  a  10  mV  electric  field  pointing  downward  (-y  di¬ 
rection),  and  this  results,  as  may  be  expected,  in  an  asymmetric  development  of  the  instability.  On 
the  right  edge  of  the  bubble,  the  background  electric  field  adds  constructively  with  the  ambipolar 
field  which  develops  to  counter  and  ion  diamagnetic  drift.  This  leads  to  a  larger  clockwise  electron 
drift,  at  least  on  the  right  hand  wall.  On  the  opposite  edge,  the  fields  add  destructively,  and  the 
electrons  experience  the  drift  of  an  effectively  much  weaker  density  gradient.  Therefore  the  insta¬ 
bility  is  nearly  stabilized  on  the  left  wall,  and  any  growth  occurring  on  the  right  edge  is  undone 
on  the  left  edge  as  the  waves  propagate  around  the  bubble.  Unfortunately  this  configuration  is 
transient  as  the  periodic  boundary  conditions  cannot  maintain  the  intermediate  scale  gradient  and 
the  stronger  instability  facilitates  faster  dissipation  of  the  initial  structure. 
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Figure  6:  Bubble  with  Efield 


3  Summary 

The  goal  for  this  project  was  to  investigate  the  small  scale  dynamics  of  CEIS.  We  performed  a 
series  of  2  and  3  dimensional  PIC  simulations  investigating  the  effect  of  strong  density  gradients 
on  a  collisionless  plasma.  The  resulting  drift  wave  instabilities  produced  irregularities  of  similar 
length  scale  as  those  observed  with  coherent  echo  radars.  At  a  maximum,  a  gradient  length  scale 
of  22  nr  was  required  to  produce  the  small  scale  irregularities,  and  all  larger  values  resulted  in 
stable  dynamics.  This  gradient  length  scale  represents  a  steeper  gradient  than  the  steepest  value 
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(approximately  50  m)  typically  observed  with  in  situ  measurements.  We  conclude  that  a  simple 
drift  wave  instability  is  sufficient  to  explain  the  development  of  small  scale  irregularities  as  they 
quickly  dissipate  the  steepened  structures  that  drive  their  growth.  Our  investigations  using  3D 
simulations  produced  dynamics  similar  to  the  2D  simulations  in  the  plane  perpendicular  to  the 
magnetic  field.  This  validates  the  results  of  2D  simulations,  at  least  for  the  early  development 
of  the  instability.  This  work  supports  the  theory  that  large  scale  CEIS  dynamics  progressively 
produces  steeper  and  steeper  structures  on  the  intermediate  scales,  which  ultimately  are  dissipated 
through  a  drift  wave  instability. 
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